How to derive and compute 1,648 diagrams 
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We present the first calculation for many-electron atoms complete through fourth order of many- 
body perturbation theory. Owing to an overwhelmingly large number of underlying diagrams, we 
developed a suite of symbolic algebra tools to automate derivation and coding. We augment all-order 
single-double excitation method with 1,648 omitted fourth-order diagrams and compute amplitudes 
of principal transitions in Na. The resulting ab initio relativistic electric-dipole amplitudes are in an 
excellent agreement with 0.05%-accurate experimental values. Analysis of previously unmanageable 
classes of diagrams provides a useful guide to a design of even more accurate, yet practical many- 
body methods. 

PACS numbers: 31.15.Md, 31.15.Dv, 31.25.-v,02.70.Wz 
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Many-body perturbation theory (MBPT) has proven 
to be a powerful tool in physics 1] and quantum chem- 
istry Although MBPT provides a systematic ap- 
proach to solving many-body quantum-mechanical prob- 
lem, the number and complexity of analytical expressions 
and thus challenges of implementation grow rapidly with 
increasing order of MBPT (see Fig^ Indeed, because 
of this complexity it has proven to be difficult to go be- 
yond the complete third order in calculations for many- 
electron atoms (see, e.g., review jj]). At the same time, 
studies of higher orders are desirable for improving ac- 
curacy of ab initio atomic-structure methods. Such an 
improved accuracy is required, for example, in interpre- 
tation of atomic parity violation j^] and unfolding cos- 
mological evolution of the fine-structure constant a |^ . 
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FIG. 1: Number of diagrams grows rapidly with the or- 
der of MBPT. Here we show number of topologically distinct 
Brueckner-Goldstone diagrams for transition amplitudes for 
univalent atoms. We assume that calculations are carried out 
in V^~^ Hartree-Fock basis to minimize the number of di- 
agrams and we do not count "folded" 6] and normalization 
diagrams. All-order single-double (SD) excitations method 
recovers all diagrams through the third order, but misses 
roughly a half of diagrams in the fourth order. These 1,576 
missed diagrams and 72 related normalization diagrams are 
explicitly computed in the present work. 



Here we report the first calculation of transition ampli- 
tudes for alkali-metal atoms complete through the fourth 



order of MBPT. We explicitly computed 1,648 topolog- 
ically distinct Brueckner-Goldstone diagrams. To over- 
come such an overwhelming complexity we developed a 
symbolic problem-solving environment that automates 
highly repetitive but error-prone derivation and coding 
of many-body diagrams. Our work illuminates a crucial 
role of symbolic tools in problems unsurmountable in the 
traditional "pencil-and-paper" approach. Indeed, third- 
order calculation ,7,,_8] was a major research project most 
likely to have required a year to accomplish. As one pro- 
gresses from the third to the present fourth order (see 
Fig. nj there is a 50-fold increase in the number of di- 
agrams. Simple scaling shows that present calculations 
require half a century to complete. With our tools deriva- 
tion and coding take just a few minutes. Similar symbolic 
tools were developed by Perger et al. p| , ho wever their 
package is so far limited to well-studiedj^, [ifl third order 
of MBPT. In contrast, we explore a wide range of new, 
previously unmanageable, classes of diagrams. 

As an example of application of our symbolic technol- 
ogy, we compute electric-dipole amplitudes of the princi- 
pal 3p3/2 — 3si/2 and 3pi/2 — 3si/2 transitions in Na. We 
augment all-order single-double excitations method 
with 1,648 diagrams so that the formalism is complete 
through the fourth order ( see Fig^ ) . The results are 
in excellent agreement with 0.05%-accurate experimen- 
tal values [l^ . Thus our computational method not only 
enables exploration of a wide range of previously unman- 
ageable classes of diagrams but also defines new level of 
accuracy in ab initio relativistic atomic many-body cal- 
culations. Atomic units |e| = ?i = me = 47r£o = 1 are 
used throughout this paper. 

Method — A practitioner of atomic MBPT typically fol- 
lows these steps: (i) derivation of diagrams, (ii) angular 
reduction, and (iii) coding and numerical evaluation. Be- 
low we highlight components of our problem-solving envi- 
ronment designed to assist a theorist in these technically- 
involved tasks. First we briefly reiterate MBPT formal- 
ism |13| for atoms with a single valence electron out- 
side a closed-shell core. For these systems a convenient 
point of departure is a single-particle basis generated in 
frozen-core (V^~^) Dirac-Hartree-Fock (DHF) approxi- 
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mation . In this approximation the number of MBPT 
diagrams is substantially reduced lEQ- The lowest-order 
valence wavefunction is simply a Slater determi- 

nant constructed from core orbitals and proper valence 
state V. The perturbation expansion is built in powers of 
residual interaction Vi defined as a difference between the 
full Coulomb interaction between the electrons and the 
DHF potential. The n*''-order correction to the valence 
wavefunction may be expressed as 
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where is a resolvent operator modified to include so- 
called "folded" diagrams 0, projection operator Q = 

1 — l^'i"'') (^'i'^^' I, and only linked diagrams are to be 
kept. From this recursion relation we may generate cor- 
rections to wave functions at any given order of pertur- 
bation theory. With such calculated corrections to wave- 
functions of two valence states w and f , 7i"^-order con- 
tributions to matrix elements of an operator Z are 
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is a normalization correction arising due 



Here Z^"-* 

to an intermediate normalization scheme employed in 
derivation of Eq. Subscript "val, conn" indicates 

that only connected diagrams involving excitations from 
valence orbitals are included in the expansion. 

Equations (Q) and ^ completely define a set of many- 
body diagrams at any given order of MBPT. In practice 
the derivations are carried out in the second quantization 
and the Wick's theorem H] is used to simplify products 
of creation and annihilation operators. Although the ap- 
plication of the Wick's theorem is straightforward, as or- 
der of MBPT increases, the sheer length of expressions 
and number of operations becomes quickly unmanage- 
able. We developed a symbolic-algebra package written 
in Mathematica [Tj] to carry out this task. The employed 
algorithm relies on decision trees and pattern match- 
ing, i.e., programming elements typical to artificial in- 
telligence applications. With the developed package we 
fully derived fourth-order corrections to matrix elements 
of univalent systems |3| • 

This is one of the fourth-order terms from Ref. ITs 
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There are 524 such contributions in the fourth order . 



Here abbreviation Sx 



; stands for £x+£y-\-- ■ -Sz, with 



being single-particle DHF energies. Further, gijik are ma- 
trix elements of electron-electron interaction in the basis 
of DHF orbitals. The quantities gijik are antisymmetric 
combinations gijik — gijik — dijki- The summation is over 
single-particle DHF states. Core orbitals are enumerated 
by letters a,b,c and complementary excited states are la- 
belled by TO, n, r. Finally matrix elements of the operator 
Z in the DHF basis are denoted as z,, . 



The summations over magnetic quantum numbers are 
usually carried out analytically. This "angular reduc- 
tion" is the next major technically-involved step. We also 
automate this task. The details are provided in Ref. 
Briefiy, the angular reduction is based on application of 
the Wigner-Eckart (WE) theorem 18] to matrix elements 
Zij and gijki- The WE theorem allows one to "peel off" 
dependence of the matrix elements on magnetic quantum 
numbers in the form of 3j-symbols and reduced matrix el- 
ements. In the particular case of fourth-order terms, such 
as Eq. lO, application of the WE theorem results in a 
product of seven 3j-symbols. To automate simplification 
of the products of 3j-symbols we employed a symbolic 
program Kentaro developed by Takada |l9j |. 

The result of angular reduction of our sample term Q 

is 
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Here the reduced quantities (i||2||j), XL^ijkl), and 
Z]^{ijkl) depend only on total angular momenta and 
principal quantum numbers of single-particle orbitals. 

As a result of angular reduction we generate analyti- 
cal expressions suitable for coding. We also automated 
the tedious coding process by developing custom parsers 
based on Perl and Mathematica. These parsers translate 
analytical expressions into FortranQO code. The result- 
ing code is very large - it is about 20,000 lines long and 
were it be programmed manually, it would have required 
several years to develop. For numerical evaluation we 
employed a B-spline library [23 • AH the fourth-order re- 
sults were computed with a sufficiently large basis of 25 
out of 30 lowest-energy [E > mc^) spline functions for 
each partial wave through /iii/2- 

At this point we have demonstrated feasibility of work- 
ing with thousands of diagrams in atomic MBPT. Now 
we apply our computational technique to high-accuracy 
calculation of transition amplitudes in Na. 

Fourth-order diagrams complementary to single-double 
excitations method. One of the mainstays of practical ap- 
plications of MBPT is an assumption of convergence of 
series in powers of the perturbing interaction. Sometimes 
the convergence is poor and then one sums certain classes 
of diagrams to "all orders" using iterative techniques. In 
fact, the most accurate many-body calculations of parity 
violation in Cs by Dzuba et al. l2ll| and Blundell et al. 
[22] are of this kind. These techniques, although sum- 
ming certain classes of MBPT diagrams to all orders, 
still do not account for an infinite number of residual di- 
agrams (see Fig.^. In Ref. 13J we proposed to augment 
a given all-order technique with some of the omitted di- 
agrams so that the formalism is complete through a cer- 
tain order of MBPT. As in that work, here we consider 
an improvement of all-order single-double (SD) excita- 
tion method employed in Ref. 22]. Here a certain level n 
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of excitations from lowest-order wavefunction refers to an 
all-order grouping of contributions in which n core and 
valence electrons are promoted to excited single-particle 
orbitals. The SD method is a simplified version of the 
coupled-cluster expansion truncated at single and dou- 
ble excitations. 

The next step in improving the SD method would be 
an inclusion of triple excitations. However, considering 
present state of available computational power, the com- 
plete incorporation of triples seems as yet impractical for 
heavy atoms. Here we investigate an alternative illus- 
trated in Fig. n ■ we compute the dominant contribu- 
tion of triples in a direct fourth-order MBPT for tran- 
sition amplitudes. We also account for contribution of 
disconnected quadruple excitations in the fourth order. 
In Ref . [13( , we separated these complementary diagrams 
into three major categories by noting that triples and 
disconnected quadruples enter the fourth order matrix 

element Z^v via (i) an indirect effect of triples and dis- 
connected quadruples on single and double excitations 
in the third-order wavefunction — we denote this class 
as Zox3] (ii) direct contribution to matrix elements la- 
belled as .^1x2! (hi) correction to normalization denoted 
as ^norm- Further these categories were broken into sub- 
classes based on the nature of triples, so that 

(zi^}) = Zix2(T„) + Zix2(r,)+ 

V / non— SD 

^ox3(5„[r„]) + Zox3p.[T.])+ (4) 
Zox3iSc[Tc]) + Zox3iD,[Tc]) + 

Zlx2{Dnl) + Zox3 

Here we distinguished between valence (T„) and core (Tc) 
triples and introduced a similar notation for singles (S) 
and doubles (D). Notation like Sv[Tc] stands for effect of 
second-order core triples (Tc) on third-order valence sin- 
gles Sy. Diagrams D^i are contributions of disconnected 
quadruples (non-linear contributions from double exci- 
tations). The reader is referred to Ref. 13] for further 
details and discussion. 




FIG. 2: Representative fourth-order diagrams involving triple 
and disconnected quadruple excitations. 

Transition amplitudes in Na. Using our problem- 
solving environment we derived the 1,648 complemen- 
tary diagrams jl3l |. carried out angular reduction [l7j| . 



and generated Fortran 90 code suitable for any univa- 
lent system. As an example we evaluate reduced electric- 
dipole matrix elements of 3si/2 — 3pi/2,3/2 transitions in 
Na (eleven electrons) "23^ . Our numerical results are pre- 
sented in TablelU Analyzing this table we see that leading 
contributions come from valence triples T^. Similar con- 
clusion can be drawn from our preliminary calculations 
for heavier Cs atom. Dominance of valence triples (T„) 
over core triples {Tc) may be explained by smaller energy 
denominators for T„ terms. Representative diagrams for 
these relatively large contributions are shown in Fig. |21 
Based on this observation we propose to fully incorporate 
valence triples into a hierarchy of coupled-cluster equa- 
tions and add a perturbative contributions of core triples. 
Such an all-order scheme would be a more accurate and 
yet practical extension of the present calculations. 

Another point we would like to discuss is a sensitiv- 
ity of our results to higher-order corrections. In Tabled 
all large contributions add up coherently, possibly indi- 
cating a good convergence pattern of MBPT. However, 
we found large, factor of 100, cancellations of terms in- 
side the ZQx^{Sy\Ty\) class. In principle higher-order 
MBPT corrections may offset a balance between can- 
celling terms and an all-order treatment is desired. For- 
tunately, the ^ox3(<S'i;[71,]) class of diagrams (combined 
with parts of Zix2{Tv)) have been taken into accoimt in 
all-order SDpT (SD + partial triples) method 22,, 24|. 
We correct our results for the difference between all- 
order [2^ and our fourth-order values for these diagrams 
(last row of Table . These all-order corrections modify 
our final values of complementary diagrams by 15%. 



Class 


Number of 


3pi/2 - 3Si/2 


3P3/2 - 3Si/2 




diagrams 










Connected triples 






•^0x3(>S'i,[T„]) 


72 


-0.8[-3] 


-1.1[- 


3] 


Zox3{Dv[Tv]) 


432 


-2.2[-3] 


-3.0[- 


3] 


Zlx2{Tv) 


504 


-0.7[-3] 


-1.0[- 


3] 


ZnoTm {Ty ) 


72 


-0.7[-3] 


-1.2[- 


3] 


Zox3{D4n]) 


144 


-0.01[-3] 


-0.01[ 


-3] 


Zox3{Sc[Tc]) 


72 


0.06[-3] 


0.09[ 


-3] 


Zlx2{Tc) 


216 


0.03[-3] 


0.04[ 


-3] 


Total triples 


1512 


-4.3[-3] 


-6.3[- 


3] 




Disconnected quadruples 






Zox3{Dnl) 


68 


1.1 h3] 


1.6[- 


3] 


•^1X2 (-Dnl) 


68 


0.2[-3] 


0.3[- 


3] 


Total quads 


136 


1.4[-3] 


2.0[- 


3] 


Total 


1648 


-2.6[-3] 


-4.3[- 


3] 


+ <5(SDpT) 




-3.3[-3] 


-4.9[- 


3] 



TABLE L Fourth-order complementary contributions to re- 
duced electric-dipole matrix elements (3pj | |-D| |3si/2) in Na. 
Last row marked "+ (5(SDpT)" is the total value corrected 
using all-order SDpT values as discussed in the text. Nota- 
tion x[y] stands for x x 10''. 



In Table ^1 we add our complementary diagrams to 
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SD matrix elements and compare with experimen- 
tal values. Several high-accuracy experiments have been 
carried out for Na, resolving an apparent disagreement 
between an earlier measurement and calculated life- 
times [see review and references therein]. In Table UTI 
we compare with the results of the two most accurate 
experiments 0, 123. The SD method [HI overestimates 
these experimental values by 2.5 a and 2.8 cr respectively 
(a is experimental uncertainty). With our fourth-order 
corrections taken into consideration the comparison sig- 
nificantly improves. The resulting ah initio matrix ele- 
ments for both 3pi/2 — 3si/2 and 3p3/2 — 3si/2 transitions 
are in an excellent agreement with 0.05%-accurate values 
from Ref. [l^ a nd differ by \.2a from less-accurate re- 
sults of Ref. yj\ . Considering this agreement it would be 
desirable to have experimental data accurate to 0.01%. 



Singles-doubles [28] 

V /non-SD 

Total 

Jones et al. 
Volz et al. [27] 



3pi/2 - 3Si/2 



3p3/2 - 3^1/2 



3.5307 
-0.0033 
3.5274 
Experiment 

3.5267(17) 
3.5246(23) 



4.9930 
-0.0049 
4.9881 

4.9875(24) 
4.9839(34) 



TABLE II: Comparison of the calculated reduced electric- 
dipole matrix element (3p-,||-D||3si/2) of principal transitions 
in Na with experimental data. 



We demonstrated that symbolic tools can replace 
multi-year detailed development efforts in atomic MBPT 
with an interactive work at a conceptual level, thus en- 
abling an exploration of ever more sophisticated tech- 
niques. As an example, we presented the first calcula- 
tions for many-electron atoms complete through fourth 
order, a task otherwise requiring half a century to com- 
plete. Although even at this level the computed transi- 
tion amplitudes for Na indicate a record-setting ah initio 
accuracy of a few 0.01%, the calculations allowed us to 
gain insights into relative importance of various contri- 
butions and to propose even more accurate yet practi- 
cal many-body method. With an all-order generaliza- 
tion 13] of the derived diagrams we plan to address 
a long-standing problem 21, 22] of improving theoret- 
ical accuracy of interpretation of parity violation in Cs 
atom [2^ . 

We would like to thank W.R. Johnson, V.A. Dzuba, 
W.F. Perger, and K. Takada for useful discussions and 
M.S. Safronova for providing detailed breakdown of 
SDpT and SD results. This work was supported in part 
by the National Science Foundation. 
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